home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Almathera Ten Pack 3: CDPD 3
/
Almathera Ten on Ten - Disc 3: CDPD3.iso
/
fish
/
701-725
/
709
/
planets
/
moon_pos.c
< prev
next >
Wrap
C/C++ Source or Header
|
1995-03-18
|
13KB
|
246 lines
/*******************************************************************************
** MOON-POS.C - Determine the position of the moon. **
** **
** Author: Jim Cobb **
** jcobb@utah.edu.cs **
** Computer Science Dept. **
** University of Utah **
** Date: Fri Mar 25 1988 **
** **
** Copyright (c) 1988, Jim Cobb, All rights reserved. **
** **
** (This copyright prevents you from selling the program - the author grants **
** anyone the right to copy and install the program on any machine it will **
** run on) **
** **
** Reference: Jean Meeus, "Astronomical Formulae for Calculators," chapter 30 **
** Meeus states that this procedure has an accuracy of 10" in the longitude **
** of the moon, 3" in its latitude and 0.2" in its parallax. **
*******************************************************************************/
#include <stdio.h>
#include <math.h>
#ifndef PI
# define PI 3.14159265358979323846 /* Should be in <math.h> */
#endif /* PI */
#ifndef RAD
# define RAD 0.01745329251994329651 /* (PI / 180 deg) */
#endif /* RAD */
/*******************************************************************************
** Trigonometric functions in degrees. **
*******************************************************************************/
#define sind(x) sin(RAD*(x))
#define cosd(x) cos(RAD*(x))
#define MAGMOON -9.9
/*******************************************************************************
** "moon_pos" takes a time argument, t, the number of Julian centuries from **
** 1900 January 0.5 ET, and the obliquity of the ecliptic, epli. **
** It adds information about the moon's position to the logfile. **
*******************************************************************************/
void moon_pos(t,epli)
double t,epli;
{
double e,e_2,
m, /* Sun's anomaly. */
m_prime, /* Moon's anomaly. */
d, /* Moon's elongation. */
f, /* Distance of Moon from its ascending node. */
omega, /* Longitude of Moon's ascending node. */
venus_term, /* Great Venus term. */
lambda, /* Ecliptical longitude of the Moon's center. */
beta; /* Ecliptical latitude of the Moon's center. */
double phase; /* Phase of the moon */
#ifdef EUTSCH
static char namestring[] = "Mond, 180.";
#else
static char namestring[] = "Moon, 180.";
#endif
double range(),poly();
int to_int();
void geo_trans();
/*******************************************************************************
** Compute the mean values for the terms. **
*******************************************************************************/
lambda = poly(270.434164,481267.8831,-0.001133,-0.0000019,t);
m = poly(358.475833, 35999.0498,-0.000150, 0.0000033,t);
m_prime = poly(296.104608,477198.8491, 0.009192, 0.0000144,t);
d = poly(350.737486,445267.1142,-0.001436,-0.0000019,t);
f = poly( 11.250889,483202.0251,-0.003211, 0.0000003,t);
omega = poly(259.183275, -1934.1420,-0.002078,-0.0000022,t);
/*******************************************************************************
** Additive terms. **
*******************************************************************************/
e_2 = sind( 51.2 + 20.2 * t );
lambda += 0.000233 * e_2;
m += -0.001778 * e_2;
m_prime += 0.000817 * e_2;
d += 0.002011 * e_2;
venus_term = 0.003964 * sind(poly(346.560,132.870,-0.0091731,0.,t));
lambda += venus_term;
m_prime += venus_term;
d += venus_term;
f += venus_term;
e_2 = sind( omega );
lambda += 0.001964 * e_2;
m_prime += 0.002541 * e_2;
d += 0.001964 * e_2;
f += -0.024691 * e_2;
f += -0.004328 * sind( omega + 275.05 - 2.30 * t );
e = poly(1.,-0.002495,-0.00000752,0.,t);
e_2 = e*e;
/*******************************************************************************
** Bring these angles within 0 to 360 degrees. **
*******************************************************************************/
m = range( m );
m_prime = range( m_prime );
d = range( d );
f = range( f );
omega = range( omega );
/*******************************************************************************
** Calculate lambda, the ecliptical longitude of the Moon's center. **
*******************************************************************************/
lambda += 6.288750 * sind( m_prime )
+ 1.274018 * sind( 2.*d - m_prime )
+ 0.658309 * sind( 2.*d )
+ 0.213616 * sind( 2.*m_prime )
- e * 0.185596 * sind( m )
- 0.114336 * sind( 2.*f )
+ 0.058793 * sind( 2.*d - 2.*m_prime )
+ e * 0.057212 * sind( 2.*d - m_prime - m )
+ 0.053320 * sind( 2.*d + m_prime )
+ e * 0.045874 * sind( 2.*d - m )
+ e * 0.041024 * sind( m_prime - m )
- 0.034718 * sind( d )
- e * 0.030465 * sind( m_prime + m )
+ 0.015326 * sind( 2.*d - 2.*f )
- 0.012528 * sind( m_prime + 2.*f )
- 0.010980 * sind( - m_prime + 2.*f )
+ 0.010674 * sind( 4.*d - m_prime )
+ 0.010034 * sind( 3.*m_prime )
+ 0.008548 * sind( 4.*d - 2.*m_prime )
- e * 0.007910 * sind( 2.*d - m_prime + m )
- e * 0.006783 * sind( 2.*d + m )
+ 0.005162 * sind( - d + m_prime )
+ e * 0.005000 * sind( d + m )
+ e * 0.004049 * sind( 2.*d + m_prime - m )
+ 0.003996 * sind( 2.*d + 2.*m_prime )
+ 0.003862 * sind( 4.*d )
+ 0.003665 * sind( 2.*d - 3.*m_prime )
+ e * 0.002695 * sind( 2.*m_prime - m )
+ 0.002602 * sind( -2.*d + m_prime - 2.*f )
+ e * 0.002396 * sind( 2.*d - 2.*m_prime - m )
- 0.002349 * sind( d + m_prime )
+ e_2 * 0.002249 * sind( 2.*d - 2.*m )
- e * 0.002125 * sind( 2.*m_prime + m )
- e_2 * 0.002079 * sind( 2.*m )
+ e_2 * 0.002059 * sind( 2.*d - m_prime - 2.*m )
- 0.001773 * sind( 2.*d + m_prime - 2.*f )
- 0.001595 * sind( 2.*d + 2.*f )
+ e * 0.001220 * sind( 4.*d - m_prime - m )
- 0.001110 * sind( 2.*m_prime + 2.*f )
+ 0.000892 * sind( -3.*d + m_prime )
- e * 0.000811 * sind( 2.*d + m_prime + m )
+ e * 0.000761 * sind( 4.*d - 2.*m_prime - m )
+ e_2 * 0.000717 * sind( m_prime - 2.*m )
+ e_2 * 0.000704 * sind( -2.*d + m_prime - 2.*m )
+ e * 0.000693 * sind( 2.*d - 2.*m_prime + m )
+ e * 0.000598 * sind( 2.*d - m - 2.*f )
+ 0.000550 * sind( 4.*d + m_prime )
+ 0.000538 * sind( 4.*m_prime )
+ e * 0.000521 * sind( 4.*d - m )
+ 0.000486 * sind( - d + 2.*m_prime );
lambda = range( lambda );
/*******************************************************************************
** Calculate beta, the ecliptical latitude of the Moon's center. **
*******************************************************************************/
beta = 5.128189 * sind( f )
+ 0.280606 * sind( m_prime + f )
+ 0.277693 * sind( m_prime - f )
+ 0.173238 * sind( 2.*d - f )
+ 0.055413 * sind( 2.*d - m_prime + f )
+ 0.046272 * sind( 2.*d - m_prime - f )
+ 0.032573 * sind( 2.*d + f )
+ 0.017198 * sind( 2.*m_prime + f )
+ 0.009267 * sind( 2.*d + m_prime - f )
+ 0.008823 * sind( 2.*m_prime - f )
+ e * 0.008247 * sind( 2.*d - m - f )
+ 0.004323 * sind( 2.*d - 2.*m_prime - f )
+ 0.004200 * sind( 2.*d + m_prime + f )
+ e * 0.003372 * sind( -2.*d - m + f )
+ e * 0.002472 * sind( 2.*d - m_prime - m + f )
+ e * 0.002222 * sind( 2.*d - m + f )
+ e * 0.002072 * sind( 2.*d - m_prime - m - f )
+ e * 0.001877 * sind( m_prime - m + f )
+ 0.001828 * sind( 4.*d - m_prime - f )
- e * 0.001803 * sind( m + f )
- 0.001750 * sind( 3.*f )
+ e * 0.001570 * sind( m_prime - m - f )
- 0.001487 * sind( d + f )
- e * 0.001481 * sind( m_prime + m + f )
+ e * 0.001417 * sind( - m_prime - m + f )
+ e * 0.001350 * sind( - m + f )
+ 0.001330 * sind( - d + f )
+ 0.001106 * sind( 3.*m_prime + f )
+ 0.001020 * sind( 4.*d - f )
+ 0.000833 * sind( 4.*d - m_prime + f )
+ 0.000781 * sind( m_prime - 3.*f )
+ 0.000670 * sind( 4.*d - 2.*m_prime + f )
+ 0.000606 * sind( 2.*d - 3.*f )
+ 0.000597 * sind( 2.*d + 2.*m_prime - f )
+ e * 0.000492 * sind( 2.*d + m_prime - m - f )
+ 0.000450 * sind( -2.*d + 2.*m_prime - f )
+ 0.000439 * sind( 3.*m_prime - f )
+ 0.000423 * sind( 2.*d + 2.*m_prime + f )
+ 0.000422 * sind( 2.*d - 3.*m_prime - f )
- e * 0.000367 * sind( 2.*d - m_prime + m + f )
- e * 0.000353 * sind( 2.*d + m + f )
+ 0.000331 * sind( 4.*d + f )
+ e * 0.000317 * sind( 2.*d + m_prime - m + f )
+ e_2 * 0.000306 * sind( 2.*d - 2.*m - f )
- 0.000283 * sind( m_prime + 3.*f );
beta *= ( 1.
- 0.0004664 * cosd( omega )
- 0.0000754 * cosd( omega + 275.05 - 2.30 * t ));
/*******************************************************************************
** Roughly calculate the phase of the moon. **
** Added by ccount **
** Sun longitude is about 36000 * T0 + 279 **
** Phase is 0 for full moon **
*******************************************************************************/
e_2 = 36000. * t + 279 - lambda;
phase = 180. - (e_2 - ((double)(to_int(e_2/360.))*360.));
if(phase < 0.)
phase += 360.;
#ifdef EUTSCH
sprintf(namestring,"Mond, %3.0f",phase);
#else
sprintf(namestring,"Moon, %3.0f",phase);
#endif
/*******************************************************************************
** Transform to right ascension and declination, and output the result. **
*******************************************************************************/
geo_trans(lambda,beta,epli,0.,MAGMOON,"PL",namestring );
}